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Introduction 

The  purpose  of  this  report  is  to  demonstrate  the  use  of  a  split-operator  (SO)  numerical 
method  on  a  problem  of  interest  to  Waterways  Experimental  Station  (WES).  Many  of  the 
environmental  problems  of  interest  to  WES  include  advective  and  dispersive  transport  and 
reaction  processes.  Mathematical  models  for  simulating  environmental  processes  that 
include  mass  transfer,  nonlinear  reactions,  multiple  species,  and  multiple  dimensions  can 
require  significant  computational  effort.  SO  methods  can  be  used  to  significantly  reduce 
the  computational  demands  of  simulating  some  of  these  problems.  We  demonstrate  a  SO 
numerical  method  on  a  model  for  sorption  and  degradation  processes.  Data  provided  by 
WES  was  used  for  the  simulations. 

Modeling  Methods 

In  this  section  we  describe  the  mathematical  models  used  to  interpret  the  batch  and 
column  data  provided  by  WES  ( Natural  Attenuation  of  Explosives  in  Soil  and  Water 
Systems  at  DoD  Sites:  Interim  Report,  1998,  Technical  Report  EL-98-,  U.S.  Army 
Engineer  Waterways  Experiment  Station,  Vicksburg,  MS).  The  processes  that  need  to  be 
modeled  in  the  batch  systems  are  transformation  reactions  and  sorption.  Advection  and 
dispersion  are  additional  processes  that  need  to  be  considered  in  the  column  systems. 

Sorption  Equilibrium  Modeling.  The  equilibrium  distribution  of  a  solute 
between  the  aqueous  and  solid  phases  is  often  assumed  to  be  a  linear  relationship, 
particularly  at  low  aqueous-phase  solute  concentrations.  Many  systems  display  nonlinear 
equilibrium  distribution  relationships,  particularly  if  aqueous-phase  solute  concentrations 
span  several  orders  of  magnitude.  In  this  report,  a  linear  model  is  used  to  describe 
sorption  equilibrium  because  it  provides  an  adequate  description  of  the  batch  equilibrium 
data  provided  by  WES.  A  linear  model  is  given  by 

&e=KCe  (1) 

where  the  superscript  e  indicates  an  equilibrium  state,  o  is  the  solid-phase  solute 
concentration,  K  is  a  partition  coefficient,  and  C  is  the  aqueous-phase  solute 
concentration. 
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Sorption  Rate  Modeling.  A  popular  model  for  describing  sorption  rates  divides 
sorption  sites  into  rapidly  and  slowly  sorbing  fractions.  Equilibrium  between  the  aqueous 
phase  and  the  solid  phase  and  is  given  by  a  combination  of  linear  models 

(2) 
(3) 

where  the  subscripts  /  and  5  indicate  rapidly  sorbing  and  slowly  sorbing  solid-phase 
fractions,  respectively. 

A  model  for  describing  transformation  reactions  and  mass  transfer  to  rapidly 
sorbing  and  slowly  sorbing  solid  phase  fractions  in  a  batch  reactor  is  given  by 
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where  t  is  time;  ka,  k y,  and  k,  are  the  reaction  rate  constants  for  the  aqueous  phase,  the 
rapidly  sorbing  solid-phase  fraction,  and  the  slowly  sorbing  solid  phase  fraction, 
respectively;  M  is  the  solid-phase  mass;  V  is  the  aqueous-phase  volume;  and  the  subscript 
sorption  is  used  to  designate  an  unspecified  expression  accounting  for  mass  transfer.  The 
rapidly  sorbing  solid-phase  fraction  is  assumed  to  always  be  in  equilibrium  with  the 
aqueous-phase  solute  concentration.  Under  this  assumption,  substitution  of  a  linear 
equilibrium  model  into  equation  5  leads  to 

dC 
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~  kfKfC 
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sorption 


Sorption  rates  for  the  slowly  sorbing  fraction  are  governed  by  first-order  mass  transfer, 
which  can  be  written  as 

dm 

a 

sorption 


dt 


=  a{mes-ms) 


(8) 
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where  a  is  a  first-order  mass  transfer  coefficient.  Substitution  of  equations  7  and  8  into 
equation  4,  and  substitution  of  equations  8  into  equation  6  leads  to  the  coupled  equations 

(i + yKf)  Tt  =  -*-c  -  y  (4°  ‘  “ .)  ■ +  kfKrc)  w 


(10) 


Initial  conditions  consistent  with  the  modeling  approach  are 


C(t  =  0) 


Kolm 
V  +  MKf 


&s(t  =  o)  =  o 


(11) 

(12) 


where  MS0]Ute  is  the  initial  solute  mass  in  the  reactor. 

An  analytical  solution  can  be  derived  for  equations  9  through  12  because  linear 
equilibrium  models  were  used.  We  developed  a  computer  program  that  solved  these 
equations  numerically  to  allow  for  the  more  general  case  of  nonlinear  equilibrium  models. 

Column  Model.  Solute  transport  in  the  soil  column  experiments  can  be 
adequately  modeled  as  a  one-dimensional  system  with  a  uniform  flow  field.  The 
governing  system  of  equations  describing  fate  and  transport  in  the  soil  column  experiments 
is 


Rr  § =  ■ D*  S-  - ~ v.  - k-c -  fW® : :  - 1 » ■ ■)  ■ +  w) 


d2C  6C 


dt 


x  dx2 


dx 


dt 


=  a(®es 


(13) 

(14) 


Rf  =  l  +  ?±Kf 
f  0  / 


(15) 


where  Rf  is  the  retardation  factor,  Dx  is  the  hydrodynamic  dispersion  coefficient,  x  is  the 
spatial  coordinate,  vx  is  the  mean  pore  velocity,  p b  is  the  bulk  density,  and  0  is  the 
porosity.  Boundary  and  initial  conditions  for  the  column  reactor  are 


C(x  =  0,  /)  =  Ca 

dCjx  =  £,*)  _  Q 
dt 


(16) 

(17) 
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C(0  <x<L,t-0)  =  0  (18) 

q  s(0  <  x  <  L,  t  -  0)  =  0  (19) 

where  CQ  is  the  influent  aqueous-phase  solute  concentration  and  L  is  the  length  of  the 
reactor. 

An  analytical  solution  can  be  derived  for  equations  13  through  19  because  linear 
equilibrium  models  were  used.  However,  we  solved  this  problem  numerically  as  a  means 
of  demonstrating  the  iterative  split-operator  (ISO)  algorithm. 

Iterative  Split-Operator  Algorithm.  The  ISO  algorithm  is  given  by:  (1)  solving 
the  transport  portion  of  the  problem  over  a  full  time  interval,  assuming  the  reaction  and 
mass  transfer  contributions  are  known;  (2)  solving  the  reaction  and  mass  transfer  portion 
of  the  problem  over  a  full  time  interval,  assuming  the  transport  contributions  are  known; 
and  (3)  iterating  over  the  first  two  steps  in  the  algorithm  until  a  convergence  criterion  is 
satisfied. 

For  simplicity  in  describing  the  algorithm,  we  restate  equations  13  and  14  in 


compact  operator  notation  as 


s* 

II 
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(24) 

The  first  step  of  the  algorithm  is  to  obtain  a  solution  over  a  time  step  t 

e  {t,  t  +  At}  to  an 

approximation  of  equation  20.  An  approximation  of  equation  20  that 

equation  21  is  given  by 

decouples  it  from 

R'^=L'+Z-' 

(25) 
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where  LrX  is  an  approximation  to  the  reaction  operator.  If  the  Crank-Nicolson  method  is 
used  with  a  central-difference  approximation  for  Lt,  then  equation  25  can  be  rewritten  as 


R 


C-c'.fcf'+fc),'  .fcj 
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1/+1  A 
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(26) 

(27) 
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where  Lt  is  a  central  difference  approximation  for  /.,,  the  subscript  i  is  a  nodal  index,  the 
superscript  l  a  temporal  index,  and  n  is  the  number  of  spatial  nodes.  A  reasonable 
estimate  of  {LrX  j  for  the  first  iteration  over  the  time  step  l  <e  {?,  t  +  At)  would  be  the  last 


estimate  of  ^Lr] )  from  the  previous  time  step  t  e  {/,  /  -  At}.  The  linear  system 


represented  by  equation  26  is  solved  to  yield  an  approximate  solution  for  the  /+ 1  time 
step. 

The  second  step  of  the  algorithm  is  to  obtain  a  set  of  solutions  over  a  time  step  t  e 
{t,  t  +  At}  to  an  approximation  of  equations  20  and  21  at  each  node  location.  An 
approximation  of  these  equations  that  decouples  them  from  the  transport  portion  of  the 
full  problem  is  given  by 


Rf  J  =Wj+(4l),  for  7  =  1,/? 

(29) 

40  ,  v 

dt  '  =(Lr2),  for  7  =  1,7? 

(30) 

where  is  an  approximation  to  the  transport  operator  at  node  i. 

(4)  can  be 

approximated  from  the  results  of  the  first  step  of  the  algorithm  as  follows 

(^1+ 1 

[L),=R>  '  a,  '  iL»\  for,  =  1’" 

(31) 

Equations  29  and  30  are  solved  over  a  time  step  /  e  {/,  t  +  At}  for  each  node.  Methods 
such  as  explicit  Euler  or  explicit  Runge-Kutta  can  be  used  for  simple  reaction  schemes. 
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More  difficult  reaction  schemes  may  require  routines  such  as  VODE,  SODEX,  or 
DASPK.  The  new  estimates  of  the  solution  at  t  +  At  are  used  to  calculate  a  new  estimate 

of  the  reaction  operator  (t,,)  for  the  first  step  of  the  algorithm.  The  update  is  calculated 
by 

tJ+1  _  . 

(L-)rRfA foTi  =  ln  (32) 

Iteration  between  the  two  steps  continues  until  convergence  is  achieved. 


Results 

We  modeled  the  sorption  and  degradation  of  trinitrotoluene  by  two  soils.  The 
experimental  results  for  the  SM  silty  sand  and  CL  lean  clay  were  used  because  these  soils 
exhibited  the  greatest  rate-limited  sorption  of  the  four  soils  investigated.  Some  of  the 
relevant  properties  of  the  two  soils  are  given  in  Table  1. 


1  Table  1  Soil  Pronerties  fdata  Drovided  bv  WES) 

%  Sand 

%  Silt 

%  Clav 

%  TOC 

SM  silty  sand  85 

7.5 

7.5 

0.020 

CL  lean  clay  65 

20 

15 

0.162 

Sorption  Equilibrium  Experiments.  Sorption  equilibrium  studies  were 
conducted  at  a  solids/liquid  ratio  of  4  grams  to  16  milliliters  over  a  24  hour  period  (data 
provided  by  WES).  Least  squares  regression  was  used  to  fit  the  linear  equilibrium  model 
to  the  TNT  sorption  data  for  the  two  soils.  The  results  are  shown  in  Figure  1.  The 
partition  coefficients  for  the  SM  silty  sand  and  CL  lean  clay  are  0.148  and  0.320  cm  /g, 
respectively. 
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Aqueous  Phase  Concentration  (mg/L) 


Figure  1 .  Equilibrium  sorption  relationships  for  TNT 
on  the  SM  silty  sand  and  CL  lean  clay. 


Sorption  Rate  Experiments.  Sorption  rate  model  parameters  can  be  estimated 
from  the  results  of  sorption  rate  experiments  conducted  in  batch  reactors.  The  form  of  the 
equilibrium  relationship  needs  to  be  specified.  This  information  is  usually  derived  from  the 
sorption  equilibrium  experiments.  However,  there  were  great  discrepancies  between  the 
equilibrium  observed  in  the  batch  equilibrium  studies  and  the  batch  rate  studies.  The  cause 
of  these  discrepancies  is  unknown.  Because  of  these  discrepancies,  it  became  necessary  to 
simultaneously  estimate  the  equilibrium  model  parameters  (i.e.,  Kf,  and  Ks)  and  rate  model 
parameters  (i.e.,  a)  from  the  batch  rate  experiments.  Degradation  rate  constants  were 
estimated  from  the  column  data  (discussed  later). 

Sorption  rate  studies  were  conducted  at  a  solids/liquid  ratio  of  90  grams  to  360 
milliliters  over  a  240  hour  period  (data  provided  by  WES).  The  results  of  the  batch 
sorption  rate  studies  and  model  fits  are  shown  in  Figures  2  and  3.  There  is  good 
agreement  between  the  data  and  model  fits,  which  is  expected  given  that  three  parameters 
were  estimated.  The  model  parameters  are  given  in  Table  2. 


Table  2. 

Model  Parameters  from  the  Batch  Rate  Exoeriments  1 

K  (cm3/g) 

/ 

A/(cm3/g)  Ks  (cm3/g) 

a  (1/hr) 

SM  silty  sand 

1.510 

0.484 

0.731  0.779 

0.0120 

CL  lean  clay 

0.606 

0.049 

0.030  0.576 

0.0132 

Uptake  of  TNT  by  the  SM  silty  sand  was  greater  than  uptake  by  the  CL  lean  clay,  which  is 
the  opposite  of  what  was  observed  in  the  equilibrium  studies.  The  SM  silty  sand  had 
roughly  equal  components  of  fast  and  slow  sorption  sites  (f  =  0.484,  where  f  =  Kf/  K). 
Almost  all  of  the  sorption  by  the  CL  lean  clay  was  rate  limited  (f=  0.049). 
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Figure  2.  Sorption  rate  study  and  model  fit  for  TNT 
on  SM  silty  sand. 
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Figure  3.  Sorption  rate  study  and  model  fit  for  TNT 
on  CL  lean  clay. 
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Sorption  Column  Experiments.  The  results  from  two  column  reactor 
experiments  were  modeled.  Dispersivity  values  are  needed  to  accurately  model  dispersion 
processes  in  column  experiments.  Dispersivity  is  usually  estimated  from  the  results  of 
conservative  tracer  studies.  These  studies  were  not  available  so  a  dispersivity  of  0.1  cm 
was  assumed.  This  value  is  typical  for  column  studies  that  we  have  conducted.  TNT  is 
known  to  degrade  in  these  systems  so  estimates  of  the  reaction  rate  constants  for  the 
aqueous  phase,  the  rapidly  sorbing  solid-phase  fraction,  and  the  slowly  sorbing  solid  phase 
fraction  were  needed.  We  assumed  that  the  three  reaction  rate  constants  were  the  same. 
The  columns  achieved  steady-state  effluent  concentrations  during  the  sorption  phase  of 
the  experiments.  The  reaction  rate  constant  ( k )  can  be  estimated  under  steady-state 
conditions  by  the  following  equation. 

-.n(C„/Q)  (33) 

where  Css  is  the  steady-state  effluent  concentration  and  th  is  the  hydraulic  retention  time. 
Initial  attempts  to  simulate  the  results  from  the  column  experiments  using  model 
parameters  estimated  from  the  batch  experiments  were  unsuccessful.  Therefore,  the 
model  parameters  (i.e.,  Kf,  Ks,  and  a)  were  estimated  from  the  column  data.  The  relevant 
information  and  estimated  model  parameters  needed  for  modeling  the  column  experiments 
are  given  in  Table  3. 
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Table  3.  Model  Parameters  from  the  Column  Experiments 


SM  siltv  sand 

CL  lean  clav 

influent  concentration  (mg/L) 

45.85 

45.85 

pore  velocity  (cm/hr) 

0.389 

0.434 

dispersion  coef.  (cm2/hr) 

0.0389 

0.0434 

bulk  density  (g/cm3) 

1.14 

1.13 

porosity 

0.459 

0.449 

decay  coef.  (1/hr) 

8.30  x  10‘5 

1.05x1  O'4 

K  (cm3/g) 

1.091 

0.493 

/ 

0.648 

0.570 

Kf{  cm3/g) 

0.707 

0.281 

Ks  (cm3/g) 

0.384 

0.212 

a  (1/hr) 

0.0269 

0.172 

The  experimental  and  modeling  results  are  shown  in  Figures  4  and  5.  The  model 
results  are  in  good  agreement  with  the  data  for  the  sorption  phase  of  both  experiments. 
The  model  does  not  accurately  simulate  the  desorption  phase  of  the  CL  lean  clay 
experiment.  The  decay  coefficients  for  the  SM  silty  sand  and  CL  lean  clay  correspond  to 
half  lives  of  348  and  276  days,  respectively.  The  mass  transfer  coefficient  for  the  CL  lean 
clay  is  larger  than  the  mass  transfer  coefficient  for  SM  silty  sand.  This  is  consistent  with 
the  lower  TNT  uptake  by  the  CL  lean  clay  (i.e.,  the  rate  of  mass  transfer  is  inversely 
related  to  sorption  uptake). 


13 


C/Co 


Figure  4.  Column  study  and  model  fit  for  TNT  on  SM  silty  sand. 
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Figure  5.  Column  study  and  model  fit  for  TNT  on  CL  lean  clay. 
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